Draft version August 7, 2012 

Preprint typeset using I4TgX style emulateapj v. 12/16/1 1 



RELATIVISTIC MAGNETIC RECONNECTION IN PAIR PLASMAS IN THREE DIMENSIONS 

Daniel Kagan Milos Milosavljevic \ and Anatoly Spitkovsky 2 

1 Department of Astronomy, University of Texas at Austin, Austin, TX 78712 
2 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544 
Draft version August 7, 2012 

ABSTRACT 

We investigate magnetic reconnection and particle acceleration in relativistic pair plasmas with three-dimensional 
particle-in-cell (PIC) simulations of a kinetic-scale current sheet in a periodic geometry. We include a guide 
field that introduces an inclination between the reconnecting field lines and explore outside-of-the-current sheet 
magnetizations that are significantly below those considered by other authors carrying out similar calculations. 
Thus our simulations probe the transitional regime in which the magnetic and plasma pressures are of the same 
order of magnitude. The tearing instability is the dominant mode in the current sheet for all guide field strengths, 
while the linear kink mode is less important even without guide field. Oblique modes seem to be suppressed 
entirely. In its nonlinear evolution, the reconnection layer develops a network of interconnected and interacting 
magnetic flux ropes. As smaller flux ropes merge into larger ones, the reconnection layer evolves toward a 
three-dimensional, disordered state in which the resulting flux rope segments contain magnetic substructure on 
plasma skin depth scales. Embedded in the flux ropes, we detect spatially and temporally intermittent sites of 
dissipation reflected in peaks in the parallel electric field. Magnetic dissipation and particle acceleration persist 
until the end of the simulations, with simulations with higher magnetization and lower guide field strength 
exhibiting greater and faster energy conversion and particle energization. At the end of our largest simulation, the 
particle energy spectrum attains a tail extending to high Lorentz factors that is best modeled with a combination 
of two additional thermal components. We confirm that the primary energization mechanism is acceleration 
by the electric field in the X-line region. The highest energy positrons (electrons) are moderately beamed 
with median angles ~ 30° -40° relative to (the opposite of) the direction of the initial current density, but we 
speculate that reconnection in more highly magnetized plasmas would give rise to stronger beaming. Lastly, 
we discuss the implications of our results for macroscopic reconnection sites, and which of our results may be 
expected to hold in systems with higher magnetizations. 

Keywords: acceleration of particles — instabilities — magnetic fields — magnetic reconnection — plasmas — 
relativity 



1. INTRODUCTION 



Magnetic reconnection (e.g., |Yamada et al.|20i0} and ref- 
erences therein) is of interest in diverse areas of astrophysics, 
yet its mechanics remains incompletely understood. A tearing 
instability is thought to be necessary to initiate reconnection in 
reversing magnetic field configurations, but many astrophys- 
ical plasmas have collisional resistivities that are insufficient 
to facilitate its g rowth. Interpretations of space plasma mea- 
surements (e.g., Chen et al. 2008; Qier oset et al.||20fl) and 
astronomical observations, however, suggest that efficient col- 
lisionless reconnection is ubiquitous. Therefore, collisionless 
effects, which operate on plasma kinetic scales, seem to be 
required to provide the dissipation necessary for effecting a 
change of magnetic topology. Magnetohydrodynamic (MHD) 
models for the dynamics of systems undergoing magn etic re- 
connection have been available for a long time (e.g., 
& Forbes 2000 and references therein). Unfortunately, 
models do not describe the underlying nature of (possibl 
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underlying nature of (possibly mul- 
tiscale) plasma organization in the reconnection layer where 
magnetic energy is being dissipated and the assumptions of 
ideal MHD do not apply. Understanding the detailed plasma 
organization on all length scales, from the likely relatively 
small, plasma kinetic scales, to the potentially much larger 
scales on which astrophysical dynamical systems "prepare" 
reconnection sites, and where ideal MHD may be valid, is 
paramount for completing the theories of a wide variety of 
astrophysical phenomena and for interpreting space plasma 



measurements and astronomical observations. 

In an effort to develop a picture of magnetic reconnection 
from first principles, recent particle-in-cell (PIC) simulations 
have examined the dynamics of reconnection layers that start 
with a current sheet as thin as the plasma skin depth. PIC 
simulations of reconnection in pair plasmas have been carried 
out in two spatial dimen sions (Zenitani & Hoshino 200l| |2007| 
IJaroschek etaLl|2004l |Bessho & B hattacharjee 2005,j2UU71 
20 1 01 120 1 2[ |Daughton &T Carimabadi 2007, Hesse & Zemtam| 
2007, Drake et al. 2010, Hoshino 2012) and three dimensions 
(jZenitani & Hoshino|2005| |2008| |Yin et al.|2008] |Liu et al | 
|2011||Sir oni & Spitkov sky|2011|>. Among these, several (|Zen-f 
itani & Hoshino 2005 1 [20081 [B essho & B hattacharj ee |2007 [ 
|Hesse & Zenitani|2007[|Hoshino|2012| ) have investigated the 
role of departure from the idealized, exactly antiparallel recon- 
nection by introducing a perpendicular "guide" field. These 
various simulations have revealed novel forms of small-scale 
plasma self-organization that are interesting in their own right, 
but that must ultimately be related to and embedded within 



the app ropriate larger astrophysical contexts (e.g., Uzdensky 
et al. 2010). While spacecraft measurements, which can be 
done in situ, can provide direct clues how to establish this 
embedding in space plasmas, in extrasolar contexts only an 
indirect relation can be established betw een the reconnection 
process and the observed e mission (e.g., |Siron i & Spitkovsky 
|20TTl|Ceruttietal.|2012b| ). 

Common features seen in many PIC simulations of magnetic 
reconnection include the formation of chains of magnetic flux 
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ropes (in three dimensions with a guide field; otherwise, the 
common terms "islands" or "plasmoids" may still be more 
appropriate), the merging of smaller flux ropes into larger ones, 
and an energization of the plasma in the reconnection layer. In 
three dimensional simulations, kink-like and oblique modes, 
as well as secondary instabilities, can impart three dimensional 
structure to the reconnection layer. 

Typically, the simulations are initialized in the so-called Har- 
ris sheet equilibrium describing a current sheet with a thickness 
similar to the plasma skin depth. The tearing instability first 
sets in on scales of the initial current sheet thickness. Its non- 
linear development produces a chain of skin-depth-scale flux 
ropes alternating with magnetic X-lines, the three dimensional 
generalization of two dimensional X-points. In X-lines, vio- 
lation of flux freezing and magnetic line rec onnection can be 
facilitated by a pressure tensor anisot ropy (fVasyIiuna s|1975| 
see also, e.g., |Hesse & Zen itani 2007] and references therein). 
Smaller flux ropes tend to merge with each other to form larger 
ones; this gives rise to magnetic organization on increasingly 
larger spatial scales. Three-dimensional PIC simulations of 
guide field reconnec tion in electron-ion p lasmas exhibit these 
same features, e.g., |Daughton et aT] ( |201 l| l found that oblique 
modes dominated over tearing modes when guide field was 
strong. Additional effects specific to plasmas with electron-ion 
mass disparity have also been identified, but are not relevant 
for the present work. 

Energization of particles in reconnectio n layers has been in- 
vestigated in a number of PIC simulations (|Zenitani & H oshino 
200T| [20071 IJaroschek et al. (20041 [Drake et al.|2006[ 120101 
Bessho & B hattachar jeel2007l|2010||2012t|Egedal et al.|2009| 
Huang etan2010[|Okaet al.|2010byLiu et al.|2011||Sironi & | 
Spifeovsky|201 l[|Egedal et al.|2012[|Hoshino|2012[[Cerutti| 



et al.| |20l2rj "Less attention has been given to particle en- 
ergization in the general case of reconnection with a guide 
field in three dimensions (Ze nitani & Ho shino 2008). In two 
dimensional simulations, particle acceleration producing a 
nonthermal energy spectrum, an apparent power-law, is of- 



ten reported. Cerutti et al. ( 2012b| l, however, instead detect 
a new ultrarelativistic thermal component energized by the 



by 



reconnection. 



agreement with analytical considerations (e.g.,|Speiser 


1965 


Larrabee et al. 2003 


Giannios|2010 Uzdensky et al. 


2011 


Cerutti et al.|2<J12ai, 


which find that particles in the vicinity 



nearly-uniform electric field as they repeatedly cross, and are 
trapped within the converging plasma flows. Other mecha- 
nisms focusing on energetic particles that have moved from 
the X-line regi on into the flanking is lands have also been sug- 
gested (e.g., Drake et al.|2006 |2010| l. In three dimensional 
simulations, evidence for a nonthermal spectrum is less solid. 
It remains poorly understood which processes limit the en- 
ergy to which particles can be accelerated in fully dynamical, 
three-dimensional reconnection layers, and how do the particle 
energy spectrum, the degree of accelerated particle beaming, 
and the temporal evolution of the accelerated population de- 
pend on the parameters of the reconnection layer. 

In this work we employ three dimensional PIC simulations 
to investigate the evolution of current sheets in relativistic 
pair plasmas undergoing magnetic reconnection. Our simula- 
tions add to the small but growing family of three dimensional 
PIC simulation of relativistic reconnection with a guide field. 
With the intention to complement existing work, we initialize 
our simulations slightly differently than it is normally done, 



not assuming the usual Harris sheet equilibrium. Also, we 
explore a parameter regime, involving magnetic to kinetic pres- 
sure ratios of the order of unity, that has thus far not received 
sufficient attention. We observe an evolution of magnetic 
field geometry that constrains the viability of models for high- 
Lundquist-number reconnection layers in which the diffusion 
region contains a hierarchy of interacting plasmoids (e.g.JShi- 



bata & Tanuma|200T Uzdensky et al.|2010] l. The simulations 
also allow us to explore the character of particle energization 
in dynamical, fully three dimensional reconnection layers. 

The paper is organized as follows. Section[2]describes our 
methodology and simulation setup, while Section [3]presents 
the results. Section|4]discusses our findings concerning devel- 
opment of kinetic instabilities in the current sheet, as well as 
our findings on particle energization, in view of the existing 
work on these topics. Finally, Section [5] reviews our main 
conclusions. 

2. DESCRIPTION OF SIMULATIONS 

2.1. The Initial Configuration 

The spatial domain is rectangular with < x < L x , < y < 
Ly, and < z < L z . The boundary conditions are periodic in 
ail directions. The initial magnetic field is the same as that in 
the Harris equilibrium 



B = Bn 



tanh 
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-tanh 
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where Ao is the half-width of the initial current sheet and k > 
is a parameter defining the strength of the uniform guide field 
perpendicular to the opposing field, which in our simulations 
is oriented in the -y direction. The current sheets are located 
at x = L x /4 and x = 3L x /4 and carry antiparallel currents. The 
current density profile that satisfies Ampere's law is 



J = 



cB 
4ttAo 



sech 



-L x /4 
An 



-sech 



-3L x /4 
An 



y. (2) 



To ensure that Ampere's law is satisfied, if we set the particle 
density to be uniform, the particles have a spatially-dependent 
drift velocity /3, = —f3 e = (3 in the -y-direction. The current J 
is related to the velocity (3 by the relation 



J = «oec(/3,-/3 e ) = 2n ecf3, 



(3) 



where e is magnitude of the unit charge carried by electron and 
ion macroparticles. This results in the drift velocity profile 



/3 = -/3o 



sech 



-L x /4 
An 



-sech 



-3L x /4 
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(4) 



where (3q = Bo/(87rnoeAo). 

In this work, we initialize the simulation outside of pres- 
sure equilibrium with a uniform initial density no and spatially 
varying drift velocity (3. We can still attempt to relate the pa- 
rameters of our simulation to those of preceding investigations. 
In setting up initial conditions for a plasma with the magnetic 
field given in Equation (QJ, there are multiple ways in which 
pressure equilibrium can be satisfied depending on the spatial 
variation of the plasma density n e + + n e -, temperature T, charge 
drift velocity (3, reversing field strength Bq, and guide field 
strength B y , It is common to assume that the initial temperature 
is uniform and equal to 7q, and only the density varies across 
the current sheet. In practice, Harris sheets are often set up 
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with a strong excess density in the current sheet for pressure 
balance. It is common to split the particle population into two 
components, one a uniform background with density n^, and 
another spatially varying with maximum density «o at the cen- 
ter of the current sheet. The latter component ensures pressure 
equilibriu m an d carries the current in the reconnection layer. 
In Section 2.4 below, we discuss how the initial configuration 



adjusts to approximate pressure equilibrium. 
One can define the magnetic to kinetic pressure ratio via 

^mag _ B 2 

~ -Pkin 87r(n e ++n <r )T' 



(5) 



where here and henceforth we express the temperature in en- 
ergy units. Then in the isothermal Harris sheet pressure equi- 
libria used in preceding investigations, the ratio simply equals 
the density contrast in the current sheet, a = no/ n b- I n our 
simulation, a is uniform outside the current sheet and equals 



16irrioT() 



(6) 



Note that a is defined not taking into account the magnetic 
pressure of the guide field. 

We do not introduce any initial perturbation to the initial 
field geometry described here. The structure that develops is 
thus seeded by numerical fluctuations. 

2.2. Parameters 

We initialize the simulation at temperature T$ = m e c 2 . All 
the particles are drawn from the relativistic Maxwellian distri- 
bution, implying that the average kinetic energy of the particles 
is ~ 2.75 m e c 2 . We are interested in the dependence of recon- 
nection mechanics and evolution of particle energy distribution 
on the dimensionless ratio of the magnetic pressure to the par- 
ticle pressure a, and the guide field amplitude parameter k. A 
nonvanishing k indicates that the magnetic field is twisted in 
the current sheet and the plasma in the sheet center is mag- 
netized. Our initial configuration allows us to freely choose 
the values of a and k in the current sheet. We run a grid of 
nine simulations, with 0.25 < a < 2 and < k < 1, and carry 
out a detailed study of the run with a = 2 and n = 0.25. The 
chosen values of a are low because we are interested in finding 
the lower limit of magnetization at which reconnection can 
produce significant particle energization. 

Because of the low growth rate, the simulations with a = 0.25 
did not develop the tearing instability over the time period of 
the simulations and thus did not enter reconnection, there- 
fore, we do not show the results of these simulations in what 
follows. The parameters of the simulations that did undergo 
reconnection are shown in Table Q] 

2.3. Simulation Method and Resolution Requirements 

We use the r elativistic particle -in-cell (PIC) plasma code 
TRISTAN-MP (Spit kovsky|2008| > to simulate the evolution of 
a reconnection configuration in a pair plasma. In a PIC simu- 
lation, electrons and positrons are represented with macropar- 
ticles. To ensure that enough macroparticles are present to 
resolve variations in the particle density, we initialize the sim- 
ulation with 4 macroparticles per species per cell. We carried 
out two dimensional test simulations with parameters match- 
ing those in the three dimensional simulations with a = 2 and 
k = 0, but with a larger number, up to 16, of macroparticles 
per cell per species. We studied various measurable quantities 



and found little variation with the macroparticle resolution, 
with the only significant discrepancy, up to 40%, seen in the 
growth rate of the kinetic tearing instability. Since in our simu- 
lations the instability is seeded by numerical fluctuations with 
different initial amplitudes in each simulation, the discrepancy 
could potentially be ascribed to an imprecision arising from the 
coarse time discretization employed in the time differencing 
in the growth rate calculation. 

To ensure that our simulations resolve the plasma skin depth, 
we set Ax = A p /8, where the skin depth is given by 



A p = 



(7) m e c 2 
Snnoe 2 



(7) 



Here, (7) is the average Lorentz factor of particles in the simu- 
lation and the additional factor of 2 in the denominator reflects 
the fact that the electrons and positrons oscillate together. We 
carried out further two dimensional test simulations with up 
to 3 times larger number of grid cells per skin depth than in 
the three dimensional simulations. Here, the most resolution- 
dependent quantity was the maximum value of the reconnected 
magnetic field B x , which varied by 25% with resolution. 

We require that the particle drift velocity in Equation (|4} giv- 
ing rise to the current density in Equation Q be nonrelativistic. 
This places a constraint on the initial current sheet width Ao- 
Combining Equations ([2| and ([3]) with the definition of a in 
Equation d6l), it can be shown that 



Ao = V2<7(r-i) 

Ap fa 



(8) 



where T— 1 = 7b/((7) tn e c 2 ') is the ratio of the particle pressure 
to the particle energy density (which includes the particle rest 
energy); in our simulations, V- 1 « 0.300 With this estimate 
of the relation between Ao and 0q, in simulations with a < 1 we 
choose A () / A p = 2, while in simulations with a = 2 we choose 
-V)/A p = 3. This implies that 0q ss 0.35, ensuring that the drift 
velocities are not relativistic. 

2.4. Readjustment to Equilibrium 

The uniform initial density and temperature and nonuniform 
magnetic field in the initial configuration defined by Equations 
(QJ, ([3]), and Q are not in force balance, but they achieve 
approximate force balance following a brief readjustment. The 
initial adjustment to equilibrium results in a reduction of the 
current sheet thickness and an accompanying compression of 
the plasma by a factor 



/.-£>.. 



(10) 



where henceforth, A denotes the current sheet half-thickness 
after the readjustment. We further define 



TcO 



Ao 
c 



A 
c 



(11) 



to denote the light crossing times of the initial and the com- 
pressed current sheet. 



The index V is related to the temperature via 



1 

r-i 



1 



*3(0~') 



(9) 



where 9 = T/(m e c 2 ) and K is the modified Bessel function of the second kind. 
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Table 1 

Simulation Parameters and Results 



Run 




V 


h " 


Ao 11 


a 


K 




\A£ B \/£ B (%) C 




Xe ner /|A£ S |(%) e 


-v f 
/max 


S1K0 


64 


40 


40 


2 


1 





1.5 


72 


3.1 


6.9 


48.4 


S1K025 


64 


40 


40 


2 


1 


0.25 


1.0 


28 


3.1 


21 


34.3 


S1K1 


64 


40 


40 


2 


1 


1 


1.0 


4.3 


1.5 


39 


34.2 


S2K0 


64 


60 


60 


3 


2 





2.3 


45 


10 


16 


52.6 


S2K025 


64 


60 


60 


3 


2 


0.25 


2.3 


44 


9.3 


14 


54.4 


S2K025L 


128 


120 


120 


3 


2 


0.25 


2.3 


27 


9.1 


28 


58.9 


S2K025P 


64 


60 


60 


3 


2 


0.25 


2.3 


44 


12 


18 


115 


S2K05 


64 


60 


60 


3 


2 


0.5 


2.3 


25 


6.4 


17 


44.9 


S2K1 


64 


60 


60 


3 


2 


1 


2.3 


6.7 


3.2 


25 


42.8 



a The length scales L v , Ly, L-, and Ao are given in units of the plasma skin depth A p . 
b / c is the ratio of the original current sheet width to the current sheet width after the readjustment phase. 
c The ratio | A£g | /£g is the fraction of magnetic energy converted into kinetic en ergy du ring the simulation. 
d K ma /K is the particle kinetic energy fraction in energized particles (see Section [3.7.l| . 

e ^ener / 1 A£"g | is the ratio of the particle kinetic energy in energized particles to converted magnetic energy (see Section [3.7.1 
f 7max is the Lorentz factor of the highest energy particle in the simulation. 



The plasma takes ss 10r c o to adjust to a state close to pres- 
sure equilibrium. The magnitude of the compression is shown 
in Table [T] With readjustment, both the density of the particles 
and the strength of the guide field increase in the center of the 
current sheet in such a way as to balance the large pressure 
gradient that the reversing field exhibits from the center to the 
surface of the current sheet. Longer-term transients are present 
in the form of acoustic oscillations. These oscillations could 
have the effect of enhancing the tearing instability, which relies 
on a velocity inflow to the current sheet. The initial readjust- 
ment does not give rise to any significant magnetic dissipation 
or particle energization. 

2.5. Unstable Modes and the Size of the Simulation 

For a simulation to capture the physics of reconnection, its 
size must be sufficient to include the fastest growing modes of 
the important instabilities of relativistic Harris current sheets 
in pair plasmas. The known instabilities include the kinetic 
relativistic tearing instability (KTI; hereafter the tearing mode), 
the drift-kink instability (DKI; hereafter the kink mode), and 
the oblique mode which is similar to t he tearing mode. To 
estima te their wavelengths, we follow Zenitani & Hoshino 
( |2007| l. The growth rate wkti of the tearing mode with wave 
number k z in a Harris current sheet of half-thickness A and 
non-relativistic drift velocity /3 is given by 



wkti = m:0/3 3/2 £ z A[i-(M) 2 ]- 



(12) 



where is a dimensionless function of the plasma temperature 
that in the limits of a cold and a relativistically hot plasma 
equals 



b(T): 



2\/2 



7r \~fpm„c'J 



-1/2 



T <C m e c , 
T > m e c 2 , 



(13) 



and 7^ = (1-/3 2 ) l l 2 . The resulting maximum growth rate 
occurs for k z \ = 1 /y/3, corresponding to a wavelength of 
~ 10.8 A. If we use the form of b(T) appropriate at ultra- 
relativistic temperatures, the growth rate for this mode is 



^KTI.n 



:0.35 



2 



(14) 



The s ituation is more complicate d for the kink mode, be- 
cause as |Zenitani & Hoshino| ( |2007] l find, the analytical max- 



sional simulations by |Zenitani & Ho shino (2005 2008) 
|Daughton et al.| ( |2011| |. The typical fastest- growing obi 



imum value of the growth rate wdki occurs for k y \ > 1, a 
wavenumber for which kinetic effects become important in 
a thin current sheet. Simulation results in their Figure 20 
indicate that the maximum growth rate occurs at k v X s» 0.7 
corresponding to a wavelength of f=s 9 A. 

In addition to these two dimensional modes, it is also neces- 
sary to resolve oblique modes with k y ,k z that combine tear- 
ing and kink compone nts; these were identified in three di men- 

and 
lique 

mode in both of these simulations had kX ~ 0.2; this corre- 
sponds to a wavelength in both directions of ~ 30 A. 

To resolve all three types of modes, in all simulations except 
for the larger size run S2K025L, we set L v = L z = 20 Ao, large 
enough to contain at least one wavelength of the tearing mode 
and two wavelengths of the kink mode. Because the initial 
adjustment leads to a significant narrowing of the current sheet, 
several wavelengths of the fastest-growing tearing and drift- 
kink modes of the narrower current sheet are included in the 
simulations, and at least one wavelength of the oblique modes 
should also be resolved. The resulting overall length scales are 
L x = 64 A p and L y = L z = 40 A p or L y = L z = 60 A p , depending 
on the simulation, as is shown in Table [T] To ensure that little 
interaction takes place between the two current sheets in the 
periodic box, we set L x = 64 A p . The current sheets in our larger 
size simulation S2K025L are as thick as in the other ones, but 
their separation and the dimensions of the box are twice as 
large, i.e., L x = 128 A p and L x = L z = 120 A p . This simulation 
was carried out on a 1024 x 960 x 960 grid and contained a 
total of w 7.5 x 10 9 particles. The smaller simulations with 
a = 1 were carried out on a 512 x 320 x 320 grid and contained 
4.2 x 10 8 particles, whereas the simulations with a = 2 were 
carried out on a 512 x 480 x 480 grid with 9.4 x 10 s particles. 



2.6. Duration of the Simulations 

In our simulations the growth times of the tearing and kink 
modes are on the same order as and smaller than the Alfven 
crossing time of the box. To capture the physics of reconnec- 
tion in its nonlinear regime and on length scales similar to the 
box size, the duration of the simulation must be larger than 
these time scales. The relativistic Alfven velocity is given by 



v'l+r/iMr-D] 



(15) 
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The Alfven crossing time of the box ta, z = L z / va is 20tco in 
the runs with L z = 60 A p . We run all of our simulations for at 
least 8000 time steps, which amounts to 150i>o in the large 
simulation S2K025L. 



3. RESULTS 

The evolution of the reconnection layer in all simulations 
exhibits the same common that have been observed in previous 
PIC simulations of magnetic reconnection in pair plasmas. The 
tearing instability grows and produces a chain of alternating 
X-lines and flux ropes. Here, we adopt "flux ropes" to denote 
a magnetic structure sometimes referred to as "plasmoid" (or 
"island" in two dimensional treatments) in which plasma is 
pinched by a helical magnetic field which can have a braided 
structure. The flux ropes can have a finite lengths, with the field 
lines opening up and extending arbitrarily far from the rope 
axis. The flux ropes can also split into sub-ropes, which opens 
the possibility of an organization of the ropes in a flux rope 
network. In all simulations except for S1K0 where the kink 
instability disrupts the current sheet, the flux ropes merge in 
quasi-hierarchical fashion until either only one flux rope is left 
in each current sheet or, in the large simulation S2K025L, the 
flux ropes have been disrupted by a transition to a disordered, 
three-dimensional state. The period of flux rope merging is 
accompanied by fast magnetic-to-kinetic energy conversion. 
This conversion produces a tail of energized particles. 

In simulations without guide field, the kink instability also 
grows, resulting in some corrugation of the current sheet. For 
a = 2 in simulation S2K0, the kink instability does not disrupt 
the quasi-hierarchical merging process. For a = 1 in simulation 
S1K0, however, the kinking disrupts the merging and brings 
the two current sheets into contact where they can interact; this 
is accompanied by rapid conversion of magnetic to thermal 
energy but the resulting spectrum does exhibit signatur es of a 
secondary energized component. In agreement with Zenitani & 
Hoshino (2008), we find that that the presence of a guide field 
suppresses the kink instability. However we find that in the 
large simulation S2K025L, a transition to three-dimensional 
evolution still occurs, likely due to a lack of large-scale phase 
coherence in the tearing instability. 

We proceed to discuss our results in detail. In Section [3TT| 
we discuss the evolution of the global current sheet and the 
formation of what will turn out to be a network of intercon- 

we 



nected flux ropes in our largest simulation. In Section 3.2 



identify site s of magnetic reconnection within the network. In 
Section |3~3"1 we carry out a Fourier decomposition of perturba- 
tions in th e lar gest simulation and discuss their growth rates. 
In Section 3.4 we analyze the time scales ass ocia ted with the 
evolution of the flux rope network. In Section[33|we measure 
the overall reconnection rate, while in Section 3~6 we analyze 
components of the nonideal electric field as reflected in the 



generalized Ohm's law. Finally, in Section [3/7] we discuss 
the overall rate of magnetic-to-thermal energy conversion as 
well as the efficiencies, mechanisms, and properties of particle 
energization in the simulations. 

3.1. Formation and Evolution of Flux Rope Network 

Figure [T] shows a time sequence of slices at y = 5 A p of the to- 
tal current density J = | J| and the plasma number density n in 
simulation S2K025L. Outside the two evolving current sheets, 
plane-parallel acoustic waves resulting from the initial pres- 
sure imbalance can be seen in the plasma density. The waves 
traverse the computational box several times in the course of 



the simulation. The early evolution of each current sheet, if 
observed at a single value of y, is that correctly described by 
the familiar flux rope (or island, plasmoid) merging paradigm. 
Larger flux ropes produced by the merging of smaller ones 
contain substructure reflected in multiple, curved, embedded 
current sheets, each with a half-thickness ~ A p . Very simi- 
lar substructure also appears in the recent t wo-dimensional 
PIC si mulations of an electron-ion plasma by Markidis et al. 
d20T2> . 

To examine three dimensional aspects of the current sheet 
evolution, in Figure [2] we show projections into the yz plane of 
the total current density J, plasma number density «, and the 
parallel electric field magnitude En in each of the two current 
sheets in simulation S2K02 5L. The flux ropes that initially 
develop from the linear tearing instability are approximately 
parallel to the y axis, which is the direction of the initial current 
flow. However, each flux rope exhibits one or more discontinu- 
ities where the z coordinate of the rope suddenly changes. We 
interpret these discontinuities as arising from a lack of tearing 
mode phase coherence on scales larger than those in causal 
contact during the initial development of the instability. The 
lack of phase coherence should occur on scales 



Ay, Az > A c 



c 



(16) 



on which locations on the current sheet are out of causal contact 
during the first e-folding of the instability; here, ui is the linear 
growth rate of a mode growing from numerical noise. This 
sets a hard upper limit on the scale on which an instability can 
grow coherently; other constraints may limit this scale further. 

The initial departure from perfect translational invariance in 
the y direction implies that the nonlinear flux rope merging will 
itself not occur coherently. In any three adjacent flux ropes, the 
middle rope can merge with one of the flanking ropes at one 
y, and with the other rope at another y; at still other values of 
y, the three ropes may remain separate for a time being. This 
creates lateral linkage between the flux ropes; already at time 
t = 37.5 t c o, all the flux ropes are mutually linked inside their 
current sheet. The linkage implies each flux rope experiences 
a strongly y-dependent magnetic tension force. This magnetic 
tension from neighboring flux ropes leads to rope tilting and 
kinking that is very distinct in origin from that arising from the 
well-known linear oblique and kink instabilities of a current 
sheet. 

Pairs of large flux ropes emerging from two generations of 
flux rope merging, but still oriented largely parallel to the y 
axis, are often connected by minor flux ropes th at are highly 
tilted toward the z axis. Dau ghton et al.| ( |2011| l have discov- 
ered similar structures in a large three-dimensional simulation 
of electron-ion reconnection with a strong guide field. It is 
clear in Figure|2]that by t ~ 75r c o, all flux rope segments have 
substantial tilts and any semblance of translational invariance 
in the direction of the original current flow is lost. Vertices 
of the flux rope network contain highly localized, intense cur- 
rent sheets and filaments. By t ~ 122r c o, the flux ropes are 
largely completely disrupted and the current sheet contains a 
disordered network of knot- and sheet-like structures. 

The force J x B peaks at the primary X-line current sheets as 
well as flux rope perimeters, but is on average much smaller in 
flux rope interiors. This suggests that the flux rope interiors are 
organized in a state of force free quasi-equilibrium constrained 
by a nonvanishing magnetic helicity (i.e., twisting, braiding of 
the field lines). 
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Figure 1. The total current density J = | J| (left two columns) and total particle density n (right two columns) in a slice y - 
column shows times / = (37.5, 47, 56, 65.5, 75)t c o and the second column shows ; = (84.5, 94, 103, 1 12.5, 122) tvq. 



5 A p . In each pair of columns, the first 



3.2. Dissipation and Reconnection in the Flux Rope Network 

To gain further insight into the structure of the flux rope net- 
work, we examine the spatial variation of Ey , the component 
of the electric field parallel to the magnetic field. The parallel 
field vanishes for the purely inductive electric field of ideal 
MHD and thus it can be obtained purely from the non-ideal, 
reconnection electric field R defined as the difference between 
the actual and the induction electric field, 



R = E+-(vxB). 

c 



(17) 



under way (e.g., Schindler et al. 1988)pThe parallel electric 
field may also be instrumental in maintaining the pressure 
anisotropy, that is responsible for the breakin g of flux freezing 
in the reconnection layer (see, e.g., Hesse et al.|201 \\ . 

The two rightmost columns of Figure [2] show that the mag- 
nitude of the parallel electric field E\\ = |E|| | strongly peaks in 
the thin, long current sheets containing primary X-lines and 
connecting magnetic flux ropes. However prior to flux rope 
disruption and the transition to a disordered state, the parallel 
field is normally very small inside magnetic flux ropes even 
when they contain embedded, substructure-level current sheets. 
An exception are flux rope segments directly undergoing dy- 
namical transformation, such as merging or stretching, where 



In high magnetic Reynolds number plasmas, the parallel field 
is present only in magnetic diffusion regions, and in particu- 
lar, in locations where a change of magnetic connectivity is 



2 On very small scales in the simulation, Em exhibits noise in the entire 
domain due to small scale electrostatic fluctuations; in our analysis, we filter 
these fluctuations on scales < 2A„. 
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Figure 2. The total current density j = |j| (left two columns), total particle density n (middle two columns), and the parallel electric field magnitude En = |E • B|/B 
projected onto the yz plane. In each pair of columns, the left column shows a projection of the sections < x < 64 A p containing one current sheet, and 
the right shows a projection of the section 64 A p < x < 128 A p containing the other current sheet. Vertically from top to bottom, the panels show times 
t = (28, 37.5, 47, 56, 65.5, 75, 84.5, 94, 103, 1 12.5) t c q. The color scale increasing from light purple to dark red is linear except for the parallel electric field, where 
the scale is logarithmic. The color scale coverage of each projected quantity was reduced from the full variation of that quantity for enhanced visual contrast. In 
particular, clipping of the color scale for the parallel electric field excludes the range in which the parallel electric field projection is dominated by small-scale 
electrostatic fluctuations in the plasma. We boxcar smoothed E and B on scales < A p prior to computing the parallel electric field, but the fluctuations still obscure 
the variation of parallel field in the current sheet at the earliest time shown, t = 28^. 
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Figure 3. A view of the current sheet at x = L x /A at times ( = 56t c o (upper 
panel) and t = 84.5 t c q (lower panel). The volume rendering (blue color) 
shows the region with current density / > 0.2cBq/(Ait\ v ) (upper panel) and 
J > 0. 12cBo/(4ttX p ) (lower panel). The curves are the magnetic field lines 
that pass through regions of high parallel electric field En and thus are actively 
undergoing reconnection. The fraction of the field line with an enhanced value 
of Eu is rendered in green. 



very isolated, spatially and temporally intermittent sites of 
significant Eu occasionally appear. The flux rope interiors are 
not undergoing pervasive, steady reconnection, in spite of their 
complex magnetic substructure. The intermittent En is also 
common after the flux rope network has become disordered, 
indicating that the disordered network remains active by de- 
veloping transient, localized dissipation events in the braided 
network of disrupted flux ropes that effect further evolution of 
magnetic connectivity. 

Reconnection intermittency is also evident in Figure[3] show- 
ing views of the current sheets t = 56 t c q, during the initial flux 
rope merging phase, and at t = 74.5 t c q, after the flux rope has 
become disordered. The field lines shown are those passing 
through the sites where En is significant, where the field line 
color is shown in green and field line geometry is that of an 
X-line. The reconnected magnetic field snaps back to become 
incorporated into that of the flux ropes. It is also noticeable 
that the flux ropes are not cylindrical^ symmetric but highly 
flattened (ribbon-like) and exhibit a clear longitudinal twist. 

A common measure of magnetic connectivity is the concept 
of magnetic helicity, but its definition on the entir e periodic 
domain presents unique challenges (e.g., Berg er|1997 and ref- 



erences therein). In spatially and temporally localized, simply- 
connected reconnection regions, a generalized helicity can be 
meaningfully defined using the Finn-Antonsen approach (e.g., 
|Berger|[l999| ). In guide-field reconnection, a change of the 
generaliz ed helicity implies a c hange of global field line con- 
nectivity (Schindler et al. 1988). Since the generalized helicity 
is locally generated by the source term -2 E • B = -2 R • B that 
is closely related to the parallel electric field, the twisting and 
braiding of the field lines interior to a flux rope is a conse- 
quence of (possibly >'-coordinate dependent) field connectivity 
change in the primary X-line region)^] 

3.3. Instability Mode Analysis 

To identify modes present in the simulation, we carried 
out Fourier decomposition of the electric and magnetic field 
for wave vectors k parallel to the plane of the current sheet, 
k x = 0, in simulation S2K025L. We search for signatures of 
the teari ng and kink instability , as well as for obl i que m odes, 
see, e.g., |Daughton et al | ( (20TT) and |Baalrud etaT] pOT2] i . The 
tearing instability has k ]| z and creates a perturbation in B x 
and E y . The kink instability has k || y and in the linear order 
creates a perturbation in E z if a guide field is present. The 
oblique mode is similar to the tearing mode, except that the 
wave vector is tilted in the yz plane, that of the initial current 
sheet. In Figure [4] we show Fourier power spectra of the 
magnetic field component B x in the linear phase at t = 19t c q 
and at the end of the fully nonlinear flux rope merging phase 
at t = 37.5 Tcq. We also show the corresponding growth rates 
evaluated over an interval of Af = 9 t c q preceding each of the 
two times, 



w Bl (k,f) 



ln|g. r (k,Q/g jt (k,r-Af)| 
Af 



(18) 



In the figure, the growth rates are multiplied by the light 
crossing time of the compressed current sheet t c w 0.43 t c q. 
The modes with fastest growth rates are clustered around the 
k z = ±2ir(L z /9)~ l tearing mode at t = 19r c o, but the fastest 
growth then shifts to longer wavelengths and oblique direc- 
tions by t = 37.5 t c( ). 

At the end of the linear phase, the power spectrum for B x 
and E y peaks at (27r) -1 (Jby,^) = (0,±(L,/9) _I ), which is the 
tearing mode, consistent with the presence of ~ 9 flux ropes in 
each current sheet. The peak, however, is broadened, likely by 
an initial lack of phase coherence on scales given by Equation 
For the tearing mode, ojkji ~ 0.083 r" 1 . From this we 



expect lack of phase coherence on scales Ay, Az 3> A CO h,KTi ~ 
5.2 Ao ~ 16A p , which can explain the observed broadening. 
However the broadening also allows the possibility that an 
authentic oblique component is present, too. At the end of 
the flux rope merging phase, which we define as the time 
t ~ (50-60)t c o when the orderly flux rope merging gives 
way to a more disordered evolution, the peak has shifted to 
longer wavelengths (L,/2) _1 < (27r) _1 |/fc z | < (L z /6)~ l with sig- 
nificant additional power occurring in the oblique direction, 
for {2ir)~ l k y < (L y /3)~ l . We emphasize, again, that the oblique 
power may be an outcome of phase decoherence, as well as 
of secondary instabilities developing in the nonlinear regime, 
rather than of a genuine primary oblique instability mode. 

Turning to our search for the kink mode, we find that E z 
develops nonvanishing albeit weak k || y power at t > 40r e o 
with the peak wavelength corresponding to (2ir)~ l k y = (L y /3) _1 . 

3 The sign of E ■ B should be uniform in each individual X-line region. 



MAGNETIC RECONNECTION IN THREE DIMENSIONS 



9 




-15 t 



1 5 t 




1 n- 2 




1 








5 E 








-* 2 








-5 E 

-10I 




io-= 












-15 P.. 




10~ 5 





5 10 

k„ 



15 



5 10 
k„ 



1 5 



' 5 



' 



-'0 



-'5 



1 5 
10 
5 

-5 
-10 
-15 



CQ 0.101 

3 



10 15 



10 15 



k„ 



Figure 4. Fourier mode amplitudes (top panels) and growth rates 
(bottom panels) of B x in one of the current sheets at times ; = 
19t c o (left panels) and t = 37.5 t c q (right panels) in simulation 
S2K025L. The amplitude was computed with the Fast Fourier Trans- 

form (FFT) via B(k) = (N,N x N z T l £~> =1 , ^ V [-2m(iyk y /N y + 

izkz/Nz)] Bfe, i y ,i z ), where N x , N y , and M are dimensions of the computational 
grid, and the initial reversing field Bo = 0.1 in these units. The components of 
the wave vector are expressed in units of 2tt divided by the computational box 
size along the relevant direction. We only show growth rates for modes with 
amplitudes \BJ)s)\ > KT 6 . 



However, most of the power in E z is still along k || z. Other 
components of E and B are devoid of power along k || y. This 
suggests that the simulation exhibits no evidence for a linear 
kink mode and that any variation with y emerges in the non- 
linear development of the tearing mode. We compare these 
results to previ ous analyses of tearing, kink, and oblique modes 
in Section l4~T1below. 



3.4. Flux Rope Merger Timescales 

Once the tearing instability has produced a collection of 
nonlinear flux ropes at a time which we will denote with tnl, 
the ropes begin to merge. The merging time should scale with 
the island separation Az as 



' merge 



Az 
XVa' 



(19) 



where x is a dimensionless coefficient encapsulating the dy- 
namics of interaction between flux ropes. Then, the number 
of flux ropes per unit length along the current sheet, which we 
denote with /i, can be estimated as 



1 



\v a (f-r NL ) 



(20) 



By visually counting the number of flux ropes in each current 
sheet in simulation S2K025Lin Figure |2j we found that 



/' : 



1.2 



<r() 



•16 



(21) 



suggesting that tnl ~ 16r c o in that simulation. Comparing 
equations (20 1 and \2\) we estimate that \ ~ 0.4. Given that 
the reconnection outflow velocities reach only about a half of 
the Alfven velocity in the simulation, it is not surprising that 
the characteristic velocity associated with flux rope merging 
\v a is only a fraction of the Alfven velocity. 

In reconnection configurations starting from a one- 
dimensional current sheet lacking any y and z dependence, 
the three-dimensionality manifested in an emerging y depen- 
dence can arise from several effects which include an intrinsic 
obliqueness of the tearing modes, the development of kink 
modes, and a lack of coherence on scales on whic h the phase of 
the tearing mode is uncorrelated (see Section |2~5) >. In reconnec- 
tion regimes in which oblique and kink modes are suppressed, 
flux rope dynamics acquires three dimensional character when 
the flux rope separation ~ pT 1 starts substantially exceeding 



the coherence length A con ,KTi- From equations ( 14 1, ( 16 1, and 
d20b, this happens at times 



t 3> T3D.coh ~ T NL + 



A 



0.35 xva/? 3 / 2 ' 



(22) 



yielding an estimate T3D. co h ~ 30r e o for simulations S2K025 
and S2K025L. Indeed, the flux rope geometry in Figure [2] 
seems to acquire three-dimensional character at t > (40 — 
50)t C (). This is consistent with the rise of power at oblique 
wave vectors k y > in Figure [4] at similar times. This analy- 
sis suggests that the reconnection layers acquire three dimen- 
sional structure on relatively short time scales regardless of the 
growth rates of the oblique and kink modes. 

3.5. The Reconnection Rate 

We now discuss the reconnection efficiency at primary X- 
lines in the early stages of the evolution of the current sheet. 
In this regime, the reconnection site can be described using 
a typical two-dimensional model of reconnection in which 
two oppositely oriented fields are separated by a thin current 
layer in which the field reverses. In this layer, magnetic field 
lines diffuse across the plasma to reconnect at one or more 
X-lines. Magnetized plasma approaches the central plane of 
the layer, toward an X-line, with an asymptotic inflow velocity 
Vi n , which is also known as the reconnection velocity. After 
passing the X-line, plasma is expelled from the vicinity of 
the X-line to either side at the outflow velocity v out , which 
is typically assumed to equal the Alfven velocity Va- The 
orientation of the reconnected magnetic field is approximately 
perpendicular to that of the original field. 

Defining 8 and L to be the thickness and length of the current 
sheet, the conservation of mass from the inflow to the outflow 
requires 

(23) 



8 
L 



Vin 



The inflow velocity can be related to the electric field in the 
reconnection region Err which is approximately uniform and 
equal to the dynamical electric field the plasma inflow just 
outside the current sheet 



Err - -(v in xBo). 



(24) 
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Because the inflow velocity and the reversing magnetic field 
are perpendicular, the dimensionless reconnection rate r rec , 
defined as the ratio of the inflow and the outflow velocity, and 
the electric field in the reconnection region, can be related via 



(25) 



Vout (v A /c)B 



This calculation of the reconnection rate makes the common 
assumption that the outflow velocity is equal to the Alfven ve- 
locity in the inflow region. In our simulations, this assumption 
is not valid as the outflow velocity is smaller than the Alfven 
velocity. 

In simulation S2K02 5, we find that indeed, the total elec- 
tric field in the reconnection region is approximately uniform 
across the current sheet. Given that the Alfven velocity is 
va ~ 0.7c in this simulation, the dimensionless reconnection 
rate calculated from the electric field following Equation ( f25] > 
is r rec 0.05 at t = 37.5 r c o, in the middle of the initial flux 
rope merging phase. This is similar to the peak reconnec- 
tion rates of about 0.07-0.1 found in other three dimensional 



PIC simulations with higher magnetizations (e.g., |Zenitani & 
Hoshino 2008; L iu et al.|201 1) >. This comparison assumes that 
reconnection at the largest X-line can be described accurately 
by a two dimensional reconnection model. The inflow and out- 
flow velocities at f = 37.5r c o are v; n s» 0.02c and v out f=s 0.27 c, 
respectively. The peak plasma outflow velocity slightly later 
in the simulation reaches v max ~ 0.38 c, just above half of the 
Alfven velocity, but since the Alfven velocity is not ultrarel- 
ativistic, it is not surprising that we are not finding ultrarela- 
tivistic flows in any of our simulations^] If the dimensionless 
reconnection rate is calculated directly from the inflow and 
outflow velocities, it is r rec « 0.08. The difference between the 
outflow velocity and the Alfven velocity, which are equal in 
idealized two dimensional reconnection models, is likely due 
to the dynamical, non-steady state nature of the plasma flow in 
our simulations. 

3.6. Generalized Ohm 's Law Analysis 

A central question in magnetic reconnection is the nature of 
the process that facilitates the violation of flux freezing in the 
diffusion region. The total electric field E in a collisionless 
plasma is related to the magnetic field B and to moments of 
the particle distribution fun ction of particle species s by the 
generalized Ohm's law (e.g., |Krall & Trivelpiece|1973 1 



E+-(v. s 
c 



XB: 



1 



VP, 



q s n s 

+ m l ( d(p s 
q s \ dt 



+ v. 



•v< Ps 



(26) 



where (v s ) and (p s ) are the average velocities and momenta 
of particles, P s is the pressure tensor, n s is the particle number 
density, and q s and m s are the particle charge and mass. Aver- 
aging over the two species s = {i,e} we recover the non-ideal, 
reconnection electric field R on the left hand side of Equation 
(26). In collisionless plasmas with long-lived magnetization 



(counterexamples being the plasmas excited by kinetic instabil- 
ities, such as the filamentation instability), the pressure and the 
inertial terms on the right hand side of Equation J26jl are small 



To obtain an ultrarelativistic flow, one mu st have 7 a = H 
(va/c) 2 ]" 1 ' 2 2> 1, but this may not be sufficient, e.g., Zenitani & Hoshino 
(2008) found maximum flow speeds of v ut & (0.7— 0.S)c even in their two- 
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Figure 5. The variation of the terms in the generalized Ohm's law of the 
electric field across the current sheet near the largest X-line at t = 47r c o. The 
center of the current sheet is at x = 0. The lines represent the measured electric 
field Ey (solid, black), the induction component (dotted, red), the pressure 
component (dashed, green), the inertial component (dot-dashed, blue), and the 
total (pink, solid) of the three calculated electric field terms. The generalized 
Ohm's law mandates that the latter should equal the actual electric field, which 
is the case apart from fluctuations arising from numerical discreteness. The 
Ohm's law components are smoothed with a Gaussian kernel of standard 
deviation 0.33A p . It is clear that at the X-line, the pressure term dominates 
over the inertial term. 



almost everywhere. Both of these terms, however, can become 
important in the diffusion region of a magnetic reconnection 
layer. 

To determine which terms give rise to R, we calculate the 
terms on both sides of Equation ( 26 1 as they vary across the 
X-line in simulation S2K02 5. We find no significant recon- 
nection electric field in the x or z directions; there is only a 
reconnection field in the y direction, consistent with two di- 
mensional models of reconnection. Meanwhile, during the 
linear growth phase of the tearing instability, spatial gradients 
are present only in the x and z directions, implying that the 
off-diagonal xy and zy components of the pressure tensor con- 
tribute to the first term on the right hand side of Equation (26 1. 
In particular, we find that most of the pressure term is proviHed 
by dPxy/dx. FigureBlshows the y component of various terms 
contributing to the electric field. The pressure term is domi- 
nant, and the total of the inductive, pressure, and inertial terms 
approximates the electric field, with departures arising from 
shot noise. 

Th ese results are consistent with previous investigations 
(e.g.JBessho & Bhattacharjee|2005|[2007||Schmitz & Grauer| 
2006) that found that the spatial variation of off-diagonal terms 
of the pressure tensor was the main contributor to the recon- 
nection electric field. The results of our Ohm's law calculation 
with a moderate guide field, combined with the similar results 



of other simulations with strong guide field (Che et al. 201 1 
nonrelativistic simulations) and no guide field (Li u et al.|2011 



dimensional simulations with high magnetization and no guide field. 



relativistic simulations), constitute evidence that for moderate 
magnetizations and all guide field magnitudes, the reconnec- 
tion electric field is produced by the pressure term in three as 
well as in two dimension s. In an exception to this agreement, 
|Hesse & Zenitani| ( |2007| l found that the inertial term became 
important in highly magnetized relativistic current sheets with 
a guide field, which is a regime we do not probe. 

3.7. Particle Energization 

We turn to the topic of magnetic to kinetic energy conver- 
sion and particle energization. Figure [6] shows the evolution 
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Figure 6. Evolution of the total magnetic energy (red lines), particle kinetic 
energy (green lines), and the energy in the electric field multiplied by 20 for 
clarity (blue lines), all normalized to the total initial energy, in runs S2K025L 
(solid lines) and S2K025 (dashed lines). The time is expressed in units of the 
light crossing time of the current sheet t c q. 

of the total magnetic field, electric field, and particle kinetic 
energies in simulations S2K02 5 and S2K02 5L. Energy was 
conserved to within 0.3% in all the simulations. The initial 
kinetic and electromagnetic energy fraction differs in the two 
simulations because the initial current sheet occupies a frac- 
tion of the volume of the simulation box that is twice as large 
in simulation S2K025 as in S2K025L. The trends, however, 
are consistent taking the volume difference into account. An 
initial reversible fluctuation of the magnetic-to-kinetic ratio 
during the first ~ 20r f o is associated with the readjustment 
to pressure equilibrium. Dissipative conversion of magnetic 
to kinetic energy seems to first occur in a couple of bursts at 
t ~ (30-40)t c o and t ~ (50 — 60)t cu . Then, the system en- 
ters a phase, lasting until the end of the simulation, in which 
dissipation is relatively steady. The two bursts seem to coin- 
cide with initial tearing instability and the formation of ~ 9 
flux ropes (magnetic islands) across the simulation box per 
current sheet, and the subsequent flux rope merging, yielding 
~ 3 flux ropes per current sheet by t ~ 60r c o, after which the 
overall geometry of the reconnection layer becomes fully three 
dimensional. By t = 150 Tco, about 40% and 20% of the initial 
magnetic energy is converted into particle kinetic energy in the 
smaller and larger simulation, respectively. 

In Table [II we show \A£b\/£b, the fraction of magnetic 
field energy that is converted to particle kinetic energy in each 
simulation. Here, we define £ B to be the initial energy in the 
total magnetic field, which includes both the reversing field 
and the guide field, and A£b is the change of the total magnetic 
energy from the beginning to the end of the simulation. The 
energy in the reversing field only is a factor of w (1 + k 2 )~ 1 
smaller than £b if we ignore the reduction of the reversing 
field in the current sheets. Independent of this correction, 
| A£ B | / £b exhibits a strong decreasing trend with increasing 
guide field strength k. It also exhibits a weaker increasing 
trend with a; note that the very high | A£b\/£b in simulation 
S1K0 is the result of destructive current sheet interaction. In 
the simulations with n > 0.5, and also in simulation S1K025 
with k = 0.25 and a = 1, the conversion of magnetic to kinetic 
energy is substantially weaker than in the other simulations. 
We also note that the rate of development and evolution of the 
flux rope network exhibits similar behavior in the sense that 
flux ropes of a given size assemble later (~ 80t cu vs. ~ 40t c o) 
in the k > 0.5 and S1K025 simulations. 
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Figure 7. Particle energy spectra in the six smaller size simulations at time 
t = 150 7^0 (top panel; see legend) and in the large size simulation S2K025L 
at five different times (middle plan; see legend) where for reference we in- 
clude the t = 150 Tea spectrum from the corresponding smaller size simulation 
S2K025 (solid line). The bottom panel compares the final spectrum to a a 
model containing three thermal populations at different temperatures (red line; 
see text), including the spectra of the three individual thermal components 
(green lines). The spectra and uncertainties are calculated from a random 
sample containing 5% of the particles in the simulation. 



The variation of | A£b \ / £b and of the overall evolution rate 
with k is likely a consequence of the decrease of plasma com- 
pressibility wi th increasing guide field strength (e.g., |Zenitani| 
|& Hesse|2008| and references therein). Since most of the mag- 
netic to particle kinetic energy conversion takes place during 
the relaxation of reconnection field lines to their final equilib- 
rium in flux ropes, the less compressible case of a strong guide 
field leads to larger flux ropes and smaller energy conversion. 

3.7.1. Energization Efficiencies 
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FigurelT] top panel, shows the energy spectrum at a reference 
time t = 150 t c o in six of the smaller size simulations. Each 
of the spectra exhibits a peak in agreement with the initial 
thermal distribution as well as a new tail resulting from particle 
energization extending to maximum Lorentz factors in the 
range 7 max ^30-50 (for the precise values at the end of the 
simulations, see Table [TJ. The simulations with the strongest 
particle energization are those with higher magnetization a = 
2 and zero or moderate guide field n < 0.25. Rather weak 
energization is seen in all simulations with strong guide field 
k = 1 and in the simulation S1K025 with weak magnetization 
and moderate guide field. We observe an intermediate level 
of energization in simulation S1K0 with weak magnetization 
and no guide field; this simulation is unique in that at the 
reference time, the two current sheets have already interacted 
with each other. It is important to note that all of these spectra 
are calculated at the same reference time; therefore, these 
differences in energization efficiency might result from the 
differing energy conversion rates discussed in the previous 
section. However, we will find below that the spectrum exhibits 
the strongest evolution during the flux rope merging phase; 
because the merging phase is complete by the reference time 
in all the simulations, it seems that the differences in the level 
of energization between the runs cannot be attributed solely to 
the differences in the rate of reconnection. 

In order to quantify the d egree of particle energiza tion in the 
simulations, following Zenitani & Hoshino (2007 ), we com- 
pute the particle kinetic energy K ena contained in the difference 
between the measured spectrum and a thermal spectrum con- 
taining the same number of particles and the same energy at 
particle energies at which the measured spectrum is in excess 
of the thermal spectrum. We then calculate the ratio of ^T e nei 
to the total particle kinetic energy K, both computed at the 
end of the simulation. It is worth remarking that we do not 
go as far as |Zenitani & Hoshino| ( |2007||2008| ) to identify the 
particles contributing to K ener as a nonthermal population; the 
following section we show this tail may be better described 
with a combination of several thermal populations. Consistent 
with the trend observed in the shape of the spectrum, we find 
that K enel /K, which is shown in Table [T| is the largest in the 
simulations with the higher magnetization a = 2 and at most 
moderate guide field k < 0.25 and is the smallest in simulation 
S1K1 with lower magnetiza tion a = 1 and str ong guide field 
K = 1. Thus, in contrast with Zeni tani & Hoshino] ( 20Q8| >, we 
find that K ena /K decreases with increasing guide field. The 
reason for this dissimilarity seems to be that the kink instability 
does not compromise the development of reconnection in any 
of our simulations other than S 1K0; even in that run, a minor 
flux rope merger occurs before the current sheet is disrupted. 
Without the kink instability, more particle energization can 
take place. 

Another possibly more interesting measure of particle ener- 
gization is the ratio of the kinetic energy in energized particles 
to the magnetic energy converted to kinetic energy. This ratio 
can be written as 
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Here, the dimensionless coefficient £ is given by 
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where Oq = Tq/ (m e c 2 ) and all the quantities entering the def- 
inition of £ are evaluated at the beginning of the simulation. 
The ratio, which ranges between 6.9% and 39%, is shown in 
Tablefl] The trend here, which an increase of K eael /\A£B\ with 
guideneld strength, is approximately the opposite of that seen 
in K ena /K. The tre nd can be understood by noting that, as we 
will find in Section 3.7.3 most of the particle energization con- 
tributing to K aKr occurs in X-line regions. Meanwhile, A£g 
measures the change in magnetic energy both during recon- 
nection in X-line regions and, more importantly, during the 
subsequent contraction of reconnected magnetic field lines into 
flux ropes. Therefore the ratio K ener /\A£ B \ should indeed in- 
crease with stronger guide field because the guide field reduces 
plasma compressibility and moderates field line contraction. 

3.7.2. Energy Spectrum Evolution and Structure 

Figure [7] middle panel, shows the evolution of the parti- 
cle energy spectrum in the large simulation S2K025L. The 
spectrum evolves little during the linear tearing phase and the 
first round of flux rope merging at t < 37.5 r e o, but then it 
quickly hardens by t = 75 r c o as the flux rope merging con- 
cludes and the reconnection layer transitions to a disordered 
state. Further hardening takes place until the end of the simu- 
lation at t = 150 T f o- For comparison, we also show the particle 
energy spectrum at this time from the smaller, equivalent sim- 
ulation S2K02 5. The two spectra are consistent at Lorentz 
factors 7 < 28 indicating convergence with increasing box 
size at relatively low energies. However, the larger simulation 
has progressively more particles at still higher energies, with 
Lorentz factors reaching 7 sa 50. 

We proceed to model the entire spectrum in simulation 
S2K02 5L as is, while keeping in mind that in a still larger 
simulation, further evolution of the spectrum at the highest 
energies is likely to be expected. We experimented with com- 
posite populations containing thermal as well as power law 
components, and found that at Lorentz factors 7 < 35, a model 
containing three thermal populations, each described by a rela- 
tivistic Maxwellian, seems to work best. The model and the 
three thermal components are shown the bottom panel of Fig- 
urelT] We fix the temperature of the first population to be equal 
to the temperature of the initial plasma, T\ = m e c 2 . A least- 
squares fit yielded temperatures T2 =2.1 m e c 2 and 73 = 3.5 m e c 2 
for the second and third population, respectively. The first com- 
ponent contains 85% of the particles and 70% of particle ki- 
netic energy, the second component contains 13% of particles 
and 24% of particle kinetic energy, and the third component 
contains only 2% of the p article s and 6% of particle kinetic en- 
ergy. We show in Section [3.7.3| that the energized populations 
and especially the particles with 7 > 30 are located close to 
the primary X-lines (within about one Larmor radius) and in 
flux ropes. 

3.7.3. Energization Sites and Mechanism 

To gain insight in the nature of the particle energization pro- 
cess, in simulation S2K025 we select a number, N lmce = 1841, 
of particles reaching the highest energies, corresponding to 
Lorentz factors 7 > 32, at the end of the flux rope merging 
phase at t = 75 t c . We trace the orbits of these particles through- 
out the merging phase over the time interval 47 t c <t < 75 t c . 
This is the period during which the particles experience coher- 
ent energization. The particles generally begin in the current 
sheet and have momenta that, on average, are aligned with the 
direction of the electric force qjE in the current sheet, with a 
median inclination of w 30° from the direction of the force. 
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Figure 8. The locus of particle acceleration in simulation S2K025 in the 
interval 47 r c < t < 70r c corresponding to the primary particle energization 
phase. The contours show the logarithm of the y-averaged energy gain per 
particle per unit volume T(x,z) for the 2000 highest energy particles in the 
simulation in this interval (see text); the lowest contour corresponds to 0.04 
times the peak value of T. The underlying color image shows the projection 
in the xz. plane of the current density J v averaged over the interval; the current 
density is scaled linearly from light purple to dark red. 

During flux rope merging, the traced particle energies in- 
crease approximately linearly. However at the instance of merg- 
ing for the two largest flux ropes, which occurs at t ~ 61 r c , 
most of the particles incur an energy increment accounting for 
~ 15% -20% of the total energy gain. Following the flux rope 
merging phase, some traced particles gain additional energy, 
but others lose energy, and both the gain and the loss could be 
considered manifestations of a thermalization process. Indeed, 
the energization /T e ner/ K does not further increase after the flux 
rope merging is complete. 

To pin down the geometric location of particle energization, 
for the traced particles, we compute the total energy gain per 
particle per unit volume averaged over y via 

dy 

x6[x-x;(t)]dt—, (29) 
L y 

where x,(f) and v,(f) are the position and velocity of traced 
particle j at time t , S is the Dirac delta, and the time integral 
covers the flux rope merging phase. In Figure [8] we compare 
the energy gain T to the corresponding absolute value of the 
>'- and f-averaged current density The highest energy gain 
is clearly associated with the largest, best defined X-line on 
the lower left. Substantial energization is also detected in the 
outflow regions flanking this X-line. The total energy gain 
in the flanking regions, seen vertically above and below the 
X-line in the figure, is larger than that in the narrow vicinity of 
the X-line. 

It is puzzling that the overall energization efficiency is 
weaker in the second current sheet shown on the right in the 




Figure 9. The trajectories of six representative particles that attain Lorentz 
factors 7 > 30 in the simulation S2K025 (left panel), and the evolution of 
the Lorentz factor as a function of coordinate z. (right panel). The dashed line 
indicates the center of the current sheet. Note that many of the particles start at 
z ~ 12.5 A p and close to the center of the current sheet; this is the location of 
the large X-line on the lower left in Figure|8] The crosses indicate the particle 
position at the beginning of the acceleration phase. 

figure. In both sheets, the energization seems to be associ- 
ated with some current density maxima but not with others. 
Energization seems to prefer long, continuous X-line regions 
with thin current sheets, perhaps because these are the sites of 
plasma inflow. 

In Figure[9]we show the orbits of six representative traced 
particles belong to the current sheet showing on the left of Fig- 
ure [8] They all start in the current sheet within ~ 8 A p from the 
X-line and have initial Lorentz factors in the range 7 ~ 16-23. 
While in the X-line region, the particles oscillate across the cur- 
rent sheet on Speiser-like orbits. The particles reach Lorentz 
factors of 7 ~ 32 (this is how they were selected; see Figure 
[9] right panel) after drifting from the X-line region into the 
neighboring islands. At the final Lorentz factor, the particles 
have Larmor radii evaluated using the reversing magnetic field 
tl = jm e c 2 j (eBo) ~ 17 A p , larger than the width of the flux 
ropes in the simulation. Because the magnetic field in flux 
ropes is only somewhat enhanced compared to the reversing 
field, the accelerated particles are not easily trapped within the 
flux ropes. 

The electric field is nearly uniform across the reconnection 
region, hence the electric force acting on a particle trapped 
in a Speiser orbit is approximately constant in time provided 
that the the inclination of the particle's velocity relative to the 
direction of the electric force is relatively small. Most of the 
traced particles happen to fulfill this condition. Then, the work 
done by the electric force on the particle is independent of the 
initial particle energy. We have also checked that, we expect 
for the Speiser orbit acceleration, it is the y-component of the 
electric field that contributes the most to the energy gain. The 
average total energy gain per particle obtained by integrating 
T over the volume of the simulation is A7 sa 12. 

Because the rate of particle acceleration in the y-direction 
is constant in time, the total energy gain is limited by the time 
the particle spends in the X-line region (see, e.g., |Cerutti et al. 
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2012a I. Particles displaced from the very center of the X-line 
experience a Lorentz force due to the reconnected magnetic 
field B x that deflects them away from the X-line and toward the 
flux rope. Once a particle has left the X-line region and entered 
a flux rope, it no longer experiences coherent acceleration. The 
lack of acceleration internal to flux ropes of the type proposed 



by |Drake et al. ( |2006| ) was attr ibuted in previous works to 
the pr esence of the guide field ( Fu et al.||2006 Huang et al. 



2010). Here, in addition the presence of the guide field, we 



note that unlike the particles traced in D rake et aL] ( |2006| >, our 
energized particles already have very large Larmor radii by 
the time they transition from the X-line region into a flux rope 
and the magnetic geometry and dynamics also seem markedly 
different. 

3.7.4. Angular Distribution of The Highest Energy Particles 

If the motion of high energy particles accelerated in a recon- 
nection region is anisotropic, then the synchrotron radiation 
emitted by these particles carries angular dependence. Ra- 
diation emitted by beams of accelerated particles can come 
in and out of view of an observer; at any given time, the 
observer detects only the radiation emitted by particles with 
momenta making angles < 7 from the line of sight. Such 
beaming could explain the temporal vari ability often seen in 
astrophysical sources, such as G RBs (e.g., Zhang &" Yan|201 1" 



McKinney & Uz densky 2012), blazars (e.g., Nalewajk o et al.| 
2011[), that may be powered by reconnection. It could also 



influence the characteristics of the gamma-ray emission from 
other candidate reco nnection-powered astrophysical sources, 
such a s magnetar s (Thompson & Duncan 1995 ; Lyutikov 2003 ; 
Parfre y~et al.|2012 1 and pulsar wind nebulae (e.g.,|Lyu barsky 
& Kir k|200lf|Slroni & Spitkovsky|2011[|C¥rutti et al.|2012a) . 
We investigate the angular distribution of the momenta of the 
highest energy particles in our simu lation, and compare with 
the results of Cerutti et al. (2012b), who recently reported a 
high degree of beaming m a two-dimensional PIC simulation 
of pair plasma reconnection initialized with a ~ 40. 

Figure [TO] shows the orientations of the momenta of all of 
the electrons and ions with 7 > 30 in simulation S2K025L. 
We select the larger simulation for this analysis to include a 
larger range of three-dimensional effects to which the parti- 
cle anisotropy may be sensitive. The median inclination of 
electrons (ions) to the y-axis (negative y-axis) is ~ 30° -40°, 
hence a half of particles in each charge species occupies a 
fraction Aft/(47r) - 0.06-0.12 of the full solid angle. This 
moderate degree of beaming is a natural consequence of X-line 
acceleration by the electric field in the reconnection region, 
since the electric field is uniform and accelerates the particles 
preferentially along the y-axis. 

The degree of bea ming, however, is much smaller than in 
|Cerutti et al.| ( |2012b| , where the highest energy particles, with 
Lorentz factors 7 > 40, were strongly beamed, occupying a 
solid angle fraction as small as AQ/(47r) ~ 0.01. The depen- 
dence of the beaming on the parameters of the reconnection 
layer can crudely be understood as follows. Assuming for 
simplicity that the plasma is ultrarelativistic and highly magne- 
tized, and that the energy of the accelerated particle increases 
by a large factor, its final degree of beaming is proportional 
to 0/(4 7r) cx (A7/7; n it)" 2 , where 7^ is the initial Lorentz fac- 
tor of an accelerated particle, and A7 is the Lorentz factor 
gain during acceleration. Using A7 cx eEAt/(m e c), where 
E ~ r lec VABo/c is the accelerating electric field (Equation 25 1 



tion[T9|, and expressing the magnetic field in terms of the mag- 
netization parameter Bo cx o-iiqTq cx (9om e c 2 ) 2 /(eA p ) 2 where 
60 = To/(m e c 2 ) (Equations|6]and|7|), we obtain 
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This result is accurate in the limit 6q ^ 1 whereas Cerutti 
et al. set up their simulation with 0q = 0.15. Ignoring this 
concern, we estimate that the factor in parentheses on the right 
hand side of Equation ( 30 ) is similar in our simulation and in 



Cerutti et al. We therefore expect that the beaming solid angle 
ft I (4tt) in their simulation should be smaller by the inverse 
ratio of the magnetizations in the two simulations, which is 
(cceratti/c r s2K025L)~ 1 = 0.05. Thus, the crude expectation is that 
the particles accelerated in Cerutti et al. should more beamed 
within a solid angle over an order of magnitude smaller than 
in our simulation, consistent with the observed beaming ratio 
between the two simulations. In conclusion, this analysis 
suggests that high magnetization and large physical size of the 
X-line region can both give rise to a high degree of beaming. 

4. DISCUSSION 

4. 1 . Tearing, Kink, and Oblique Modes 

A number of studies have examined the three-dimensional 
evolution of a current sheet without a g uide field; here, we 
focus o n results relevant to pair plasmas. |Zenitani & Ho shino 
(2008 ) found that the kink mode dominated the tearing mode 
in antiparallel reconnection. In contrast, Liu et al.| (201 1 ) found 
that the initial evolution was dominated by the tearing mode, 
leading to an initial merging phase similar to that in two dimen- 
sions. After the completion of the merging phase, however, 
they found that a secondary kink instability set in and rendered 
the reconnection region turbulent and filled with plasmoids 
covering a range of dimensions. The reduced influence of 



the kink mode in the simulations of Liu et al. (201 1 1 may be 
explained by initially relativistic drift velocities of 0.82c. Ana- 
lytical calculations suggest that while the kink mode dominates 
at low drift velocities, the tearing mode be comes dominant 
for drift velocities exceeding 0.6 c ( jZenitani & Hoshino|2007| ). 
However, the growth rate of the kink mode, unlike the tearing 
mode, depends strong ly on the initial structure of the current 
sheet (Daughto n]l999| ). 

In guide field reconnection, previous three dimensional stud- 
ies have found that the kink mode is stabilized, but, u nder 
certain conditions, obl ique modes dominate the evolution. [Zen] 
|itani & Hos hino ( 2008 ) found that with a guide field present, 
the tearing mode and a sausage-like mode are combined into 



an oblique "relativistic-drift-sausage-tearing" mode. Daughton 
et al. (201 1) found that in electron-ion plasmas, oblique modes 



and Af ~ r merge ~ Az/v A is the duration of acceleration (Equa 



give rise to a network of interconnected flux ropes. In both 
cases, the nonlinear development of the oblique modes led to 
turbulence in the current sheet. In no nrelativistic electr on-ion 
simulations with a strong guide field, |Che et aL ( 201 1\ found 
that at low temperatures, the current sheet filamented and that 
the turbulent transport played a key role in the diffusion of the 
field. At higher, albeit still nonrelativistic temperatures, their 
filamentation instability became weaker and the current sheet 
did not become turbulent. 

One must be aware of the possibility that certain instabilities 
identified in simulations starting with Harris sheet equilibria 
are idiosyncrasies of the initial configuration and have little to 
do with astrophysical magnetic reconnection. Our simulation, 





Figure 10. Angular distribution of the momenta of particles with Lorentz factors 7 > 30 at times t = 56tvo (top panels), t = 66r c o (middle panels), and t = 84r c o 
(bottom panels) in simulation S2K025L shown in the Aitoff projection. The panels of the left and right show the particles near the current sheets at x = L v /4 and 
x = 3L x /4, respectively. The plotted particles were selected from a random sample containing 5% of the particles in the simulation. 



starting from a non-Harris like configuration lacking a density 
enhancement in the initial current sheet, should help discern 
artifacts of the initial configurations from genuine properties 
of reconnection dynamics. 

We have found that in our simulations in which we consid- 
ered guide fields at most equal to the reversing field, the tearing 
mode dominates the evolution at all guide field strengths. To 
understand the difference with the studies that detected a dom- 
inant kink mode, we carried out two dimensional simulations 
with a = 2 and n = in both the xy and xz planes. We found that 
the kink mode grows faster in a Harris current sheet, while the 
tearing mode grows faster starting with our initial conditions. 
Our initial configuration seems to enhance the growth rate 
of the tearing mode and inhibit the growth of the kink mode. 
The plasma inflow occurring during readjustment to equilib- 



rium may enhance growth of the tearing mode above the value 
expected in pressure equilibrium. Furthermore, our initial con- 
figuration may be endowed with a somewhat stronger drift 
velocity shear than the Harris sheet and this could stabilize the 
kink mode ( [Volponi et al . 2000 ) w hile having littl e effect on 
the tearing mode ( Roy tersht eyn & Daughton|2008| l. 

We proceed to calculate the theoretical and actual tearing 
mode growth rates in each simulation. In Table [2] we show the 
two rates which we label cjKTi.t and wkti.s, along with the initial 
particle drift velocity /3 , the initial current sheet compression 
factor / c , and the number of flux ropes per current sheet deter- 
mined by visual inspection A^fr- Because the fastest-growing 
tearing wavelength computing from the theoretical result in 
Equation ( 12 1 is ~ 10. 8 A and the box size in all by simulation 
S2K025L is 20 Aq, the number of flux ropes per current sheet 
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Table 2 

Tearing Instability Parameters 



Run 


Pa 


fc* 


N m b 


TcO^KTI,t C 


T cO^KTI,s C 


S1K0 


0.365 


1.5 


3 


0.12 


0.06 


S1K025 


0.365 


1.0 


2 


0.08 


0.04 


S1K1 


0.365 


1.0 


2 


0.08 


0.03 


S2K0 


0.344 


2.3 


4 


0.16 


0.14 


S2K025 


0.344 


2.3 


4 


0.16 


0.13 


S2K025L 


0.344 


2.3 


9 


0.16 


0.11 


S2K1 


0.344 


2.3 


2 


0.09 d 


0.07 



a / c is the ratio of the original current sheet width to the current 
sheet width after the readjustment phase. 

b Wfr is the number of flux ropes per current sheet initially formed 
in the simulation. 

c The two values of t c oUkti are me normalized theoretical (sub- 
script "t") and simulated (subscript "s") growth rates of the tearing 
mode. 

d This growth rate is based on the observed mode wavelength, 
which is significantly longer than the fastest-growing wavelength. 



produced by the tearing mode should be A^r = l-9/ c ; the num- 
ber should be twice as large in S2K02 5L. The actual number 
of flux ropes formed in the simulation is agreement with this 
prediction, with the exception of the simulation S2K1 with 
the strongest guide field. This run produced only 2 flux ropes 
pe r current sheet even th ough one expects 4. As discussed 
in |Daughton et al.| ( |20Tl) 



the strong guide field may prefer 
oblique rather than pure tearing modes, but the former grow 
the fastest at long wavelengths ~ 33 A p , about half the box size 
of the simulation, and the finite box size interferes with ability 
of arbitrary oblique modes to grow. 
The values of cjoi.t for the various runs are shown in Table 

tit should be noted that these are upper limits to the possi- 
e tearing mode growth rates, because any measured growth 
rates are found from estimates in the nonlinear regime. We 
also compute the simulated growth rates Wkti.s of the tearing 
mode in our simulations by calculating the growth rate of the 
average reconnected magnetic field perturbation B x in the cur- 
rent sheet. Growth rates of the KTI calculated with the latter 
method should be higher than those obt ained by examining a 
single Fourier mode, as we did in Section [3~3| because multiple 
Fourier modes contribute to the amplitude of the perturbation. 
We find that the theoretical and simulated growth rates were 
relatively similar in all simulations, wit h theoretical rates be- 
ing at most double the simulated rates. Zenitani & Hoshino 
( |2007[ ) identified a similar discrepancy between the theoretical 
and simulated growth rates. The growth rates also show that 
decreased guide field k, and especially higher magnetization 
er, lead to higher tearing mode growth rates both in theoretical 
calculations and in the simulations. This reinforces the effects 
of the variation in A^r, which has a similar dependence on a 
and k. 

We have shown that the tearing mode is responsible for 
producing magnetic flux ropes in most of our simulations. 
The tearing mode can also explain many of the differences 
between the simulations. The runs with higher a and smaller 
k have larger values of f c , which leads to higher tearing mode 
growth rates and the formation of a larger number of flux ropes 
A^fr- The faster evolution of runs with higher a and smaller n 
can be explained by the faster growth of the tearing mode in 
such runs. Because runs with higher a and smaller k are also 
associated with a larger A^fr, more flux rope merging can take 
place in such runs. Flux rope merging is strongly associated 
with energy transfer and particle acceleration (as we show in 



the following sections), so the larger energy transfer in runs 
with higher a and smaller n can be explained by the larger A^fr 
in those runs. This suggests that a large portion of the variation 
between runs is a result of the effect of a and k on the tearing 
mode. 

4.2. Particle Energization 

Various particle energization channels have been iden tified 
in PIC simulations of plasma reconnection (see, e.g., 



Oka 



|et al.|2010b] and references therein). Three specific loci where 
particle energization was detected include: inside or near the 
diffusion region containing reconnection X-lines, in the mag- 
netic islands (or plasmoids, flux ropes) where the reconnected 
magnetic flux accumulates, and between an X-line and the 
edge of a flanking island, where the plasma flowing out of the 
X-line first encounters a strong magnetic field gradient. 

A number of studies point to the conclusion that significant 
energization occurs near the primary (i.e., largest, or highest 
rank) X-lines. There, the electric field, which is perpendicular 
to the reversing field and is aligned with the current density 
vector in the middle of the current sheet, can accelerate parti- 
cles oscillating w i thin or across the current sheet (e.g.,|Zenitani 
| & HosMno|2^[2^|2TO||Lyubarsky & Liverts 2008 fjUzl 
|densky et al.|201l[ Bessho & Bhattacharjee |2012||Cerutti et al.| 
201 2a|bp . A variant of this mechanism involves the trapping 
of particles in secondary magneti c islands appeari ng within 
the diffusion region of an X-line ( Oka et al.|2010a i; we will 
not discuss this as we do not observe secondary flux rope 
formation in our simulations. 

Particles can also be accelerated outside of the X-lines, such 
as in primary magnetic islands. If an island contracts in the 
course of its relaxation to MHD equilibrium, the particles 
trapp ed inside it can be accelerated by a Ferm i-type process 
(e.g., |Drake et al.|2006||2l)10||Kowal et al.|201 l| i. Yet another 
location for particle acceleration is in the pileup region be- 
tween the X-line and a flanking island, which is where the 
reconnected magnetic flux accumulates. There, relativistic 
Speiser motion can combine with curvature drift along the 
magnetic field gradient to create significant acceleration (e.g., 
Hoshino et al. 2001 ; Jaroschek et al. 2004; Zenitani & Hoshino 
2007; ]Pntchett|2008||Huang et al.|2010||Liu et al.|201 1| |. 

These mechanisms can operate in reconnection sites where 
the plasma configuration, involving inflow into the diffusion 
region and outflow toward the flanking islands, is relatively 
stationary, or inside a single, autonomously evolving magnetic 
island. There may be other mechanisms that operate in course 
of spatial rearrangement and merging of magnetic islands. 
Converging islands can give rise to Fermi-type acceleration 



as particles boun ce between them (e.g., Oka et al. 2010b 
Tanak a~t al.|201 1 1. A dynamically active reconnection region 



containing many interacting islands can also allow Fermi-type 
stocha stic particle acceleration ( |Drake et al.|2010"t |Hoshino| 



|2012 1. It is worth noting that even if particle acceleration is 
not directly driven by the merging, it is normally most active 
during isla nd merging episodes in the cou rse of a reconnection 
event dJaroschek et al.|2004[ |Pritchett|200"8] l. 

The presence of a guide field has a significa nt effect on par - 
ticle acceleration in reconnection (Zenitani & Hoshino|2008| l. 
Without a guide field, the direct X-line acceleration was typ- 
ically found to be less efficient than the other mechanisms 
discussed. However, with a guide field , the effective ness of 
acceleration in the pile up regi on (Huan get al.|2010] l and in- 
side magnetic islands ( |Fu et al.|2006| l was diminished, thus 
leaving X-line acceleration of particles on Speiser orbits as 
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the dominant mechanism. Our results are consistent with this 
conclusion. 

To permit an interpretation of the nonthermal radiation spec- 
tra in observed high-energy astrophysical sources in terms of 
synchrotron and inverse-Compton radiation, a relatively hard 
power law particle energy spectrum, which places a significant 
portion of the total energy in high-energy particles, is required. 
In most PIC simulations including the present work, reconnec- 
tion produced energized particle populations, but whether the 
populations represented genuine power laws tails has remained 
a matter of interpretation. Alternatively, an energized popula- 
tion can be interpreted as a composite of one or more thermal 
sub-populations, each at a different temperature, for example. 
Two dimensional simulations typically provide the dynamic 
range to make a tentative distinction between a thermal or 
a power law spectrum in an energized population. In most 
three-dimensional simulations, however, such a determination 
is dubious. In what follows, we discuss the characteristics 
of the particle energy spectra in the literature and provide a 
comparison with our results. Then, we briefly reflect on the 
expected nature of particle energy spectra produced by systems 
experiencing magnetic reconnection. 

Most investigations involving PIC simulations of plasma 
reconnection have interpreted a section of the high energy tail 
of the particle energy spectrum as a power law dN/dln-f oc 7°. 
In the X-line region of two dimensional simulations, spectra 
with power law indices as hard as a sa -1 have been reported 
l|Zenitani & Hoshino|200T1[2007l|Jaroschek et al.|2004[|Besshol 
|& Bha ttacharjee 2007; Lyubarsky & Liverts 2008). The spec- 
trum in the whole simulation box also contains a power law 
component, but with a softer index of a ~ -2.5. In their two- 
and three- dimensional simulations o f shock-induced recon- 
nection, Sironi & Spitkovsky (201 1 ) find that as long as the 



region containing a reversing magnetic field that can undergo 
reconnection is reasonably large compared to the thickness of 
the reconnection layer, the particle energy spectrum at high 
energies is a power law with an index a = -1 .5 over a decade in 
energy. Other three dimensional investigations have identified 
relatively soft power law spectra with indices a ~ -3 both in 
the X-line region and in the entire sim ulation box ( |Jaroschek| 
|et al.|2004|rZenitani & Hoshino|2008) . 

Other investigations have interpreted particle energy spectra 
in terms of multi ple therm al and other exponentially truncated 
populations (e.g., |Oka et al.|2010b| >. In their two-dimensional 
simulation of reconnection in an initially nonrelativi stic plasma 
at temperature T = 0.15ra e c 2 , Cerutti et al. ( 20 12b} find a new 
thermal population at temperature T ~ 4m e c z . In a three- 
dimensional simulation beginning with a relativistic plasma 
at temperature T = 
population with T 
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2.3 m e c 2 is produced. In a recent work 
presen ting two-dimensional simulations without a guide field, 
Bessho & Bhattacharjee (2012 ) detect a spectrum of the form 
dN/dj oc 7" 1/4 exp(-a7 1/2 ), where a is a constant of the order 
of unity. 

Among the cited descriptions of particle populations ener- 
gized by reconnection, the spectra observed in our simulations 
bear resemblance with those involving multiple thermal sub- 
populations. For example, the spectrum at the end of simula- 
tion S2K02 5L can be modeled with three the rmal component s 
at temperatures similar to those found in Liu et al. ( 2011[ ); 
specifically, the two energized sub-populations have tempera- 
tures T = 2.1 m e c 2 and T = 3.5 m e c 2 . Neither a power law nor 



(2011 



find that a thermal 



good fit to the spectrum in this simulation. It is important to 
note that the continued evolution of the particle energy spec- 
trum we observe after the current sheet has become disordered, 
and the likely additional evolution that would be taking place 
in am even larger simulation, mean that the spectral form has 
not converged. For example, it is possible that in addition to 
the two thermal components that we have detected, additional 
such populations at still higher temperatures would appear in 
larger simulations, and that the combination of those would 
constitute a hard spectral tail, e.g., a power law. 

In PIC simulations, distinguishing between the various forms 
of energized particle spectra is difficult due to dynamic range 
limitations. Another complication is how the spectrum of the 
non-energized, yet possibly adiabatically heated background 
plasma is to be subtracted to isolate the genuine nonthermal 
component. Insight can separately be gained by examining 
test particle trajectories in magnetic field geometries model- 
ing reconnection layers. Particle ballistics in X-line as well 
as magnetic island geometries introduced at various levels of 
approximation has been investigated in numerous s t udies (e.g., 
I Zenitani & Hoshino|[200T] |La rrabee et al. 2003 ; |Bessho &| 
|Bhattacharjee|2012||Cerutti et al.|2012a i. Often, the models 
entail transport terms describing particle acceleration while 
confined in a single X-line region or an island, and other terms 
describing particle escape and the termination of acceleration. 
Then, the terminal energy spectrum of the escaping particles 
is obtained by taking the input energy spectrum and deter- 
ministically transforming it by the transport terms. The crud- 
est such models have suggested that the terminal spectrum 
could be a power law, e.g., of the form dN/dj oc 7" with 
a = -(2/7r)Z? rec /£ , RR ~ — (2/7t)c/va, where B rec is the recon- 
nected magnetic field and £rr is the elec tric field in the recon- 
nection region discussed in Section 3.5 (Zenitani & Hoshino 



the Bessho & Bhattacharjee ( 2012 1 spectral form present a 



|2001| l. Models taking a more detailed accounting of the trans 
port of particle phase space coordinates and the kinematics of 
escape, however, ca n instead imply distinctly non-p ower-law, 
softer spectra (e.g., Bessho & Bhattacharjee 2012 1. It is also 
important to note that the breaking of translational invariance 
in the direction of the cur rent, e.g., by three dimensional phase 
decoherence (see Section [3~T| ), may limit particle acceleration 
and prevent the production of the same energized particle tail 
produced in two dimensional geometries. 

Power-law spectra are generically expected in acceleration 
processes entailing stochasticity. A prime example is the lin- 
ear diffusive shock acceleration (DSA) in which the random 
outcome of particle scattering in the shock downstream deter- 
mines whether the particle will return to the shoc k upstream 
and be subjected to anothe r acceleration cycle ( |Bell||1978| 
|Blandford & Ostriker|1978] l. In contrast, the acceleration in 
an idealized time-independent and two-dimensional X-line 
region that lacks substructure in the form of secondary and 
embedded islands, should b e deterministic. Therefore the mod- 
els of Zenitani & Hoshino (2001 ), Larrabee et al. (2003 ), and 
Bessho & Bhattacharjee ( 2012) 1, can be thought of as producing 
a power-law-like spectrum only "by coincidence." Realistic 
reconnection regions should plausibly allow particles to be 
accelerated in stochastic fashion as they contain time depen- 
dence and structure on multiple scales; they may also have 
a mechanism for returning particles that have es caped an ac- 
celeration site into another such site (see, e.g ., |Drake et"aL] 
2006 20i0l |Kowal et al.|2011|poshinoj2012| i. Since in this 
work we have found that the reconnection layer transitions 
into a disordered network of interacting flux ropes, it will be 
particularly interesting to investigate, in subsequent study, if 
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this network allows for stochastic acceleration of particles in 
intermittent, secondary reconnection sites that may appear in 
the course of the flux rope network evolution. 

5. CONCLUSIONS 

In this paper, we carried out three dimensional PIC simu- 
lations of magnetic reconnection in a relativistic pair plasma 
with varying guide field strength. Plasma magnetizations, ex- 
pressed in terms of the magnetic to kinetic pressure ratio, were 
of the order of unity. The initial conditions differed from the 
usual Harris sheet configuration by not having a large den- 
sity contrast between the center of the current sheet and the 
background plasma. We investigated the growth of unstable 
kinetic modes in the current sheet, as well as the nonlinear 
development of a three dimensional flux rope network. We 
also investigated the character and efficiency of particle ener- 
gization. Our main results can be summarized as follows: 

The current sheets in all simulations develop significant mag- 
netic reconnection accompanied with conversion of magnetic 
to particle kinetic energy. With the aid of Fourier decomposi- 
tion, we ascertained that the linear tearing mode is dominant 
in the early evolution of the current sheet. We find that no 
significant growth occurs in the linear kink and oblique modes. 
The nonlinear development of the tearing mode produces a 
chain of flux ropes separated by primary X-lines. The flux 
ropes merge in hierarchical fashion whereby the merging time 
scale is proportional to the flux rope separation. During this 
phase magnetic reconnection takes place at the X-lines. We 
find that the dimensionless reconnection rates ~ (0.05-0.08) 
and the maximum outflow speeds w 0.4c ~ va/2 in our simu- 
lations are similar to those detected in other three dimensional 
simulations of reconnection in pair plasmas. We also find that 
spatial variation of an off-diagonal component of the pressure 
tensor is responsible for the breaking of flux freezing at the 
X-lines, consistent with existing results. 

While the hierarchical flux rope merging process initially 
appears similar to that found in two dimensional simulations, 
in fact it is three-dimensional from the outset. This is because 
a lack of initial phase coherence in the linear tearing mode on 
scales larger than those allowed by causality breaks transla- 
tional invariance in the direction of the initial current flow. The 
flux ropes form a topologically interconnected, dynamically 
evolving network. Dynamical interaction between neighbor- 
ing flux ropes is provided by magnetic tension forces. With 
time, the flux ropes break up into segments with more isotropic 
orientations. The strongly three-dimensional character of the 
reconnection layer seems to suggest that global reconnection 
models invoking quasi-two-dimensional plasmoid hierarchies 



(e.g.JShibata & TanumaJ2001||Fermo et al.|2010"{|Uzdensky 
|et al.|2010p require revision to account for the inter-plasmoid 
magnetic linkage and isotropization of plasmoid orientations. 

The larger flux ropes produced during flux rope merging 
contain substructure down to plasma skin depth scales which 
is reflected in embedded, twisted and braided current filaments 
and sheets. Overall, this substructure is force-free and evolves 
relatively slowly. However, isolated sites within the evolved 
flux rope network contain spatially and temporally intermit- 
tent sites characterized by strong nonideal conditions E • B ^ 
where a change of magnetic connectivity continues to take 
place even after flux rope merging has saturated on length 
scales equal to the size of the computational box. This inter- 
mittency may produce the observed variability of nonthermal 
emission in systems in which the emitting particles are ener- 
gized by magnetic reconnection. 



During the early, ordered flux rope merging phase, particles 
are accelerated to high Lorentz factors by the electric field 
in primary X-lines; the trajectories of these particles are well 
described by Speiser orbits. Particles continue to be energized 
in the later, disordered phase we identify in our largest simula- 
tion, but we leave the analysis of energization in the disordered 
regime to a subsequent investigation. 

Simulations with higher magnetization and lower guide field 
strength exhibit greater and faster energy conversion and parti- 
cle energization. The efficiency of particle energization mea- 
sured in terms of the energy in the accelerated particles per unit 
magnetic energy dissipated in the simulation is an increasing 
function of the guide field strength, which can be interpreted 
as resulting from a decreasing plasma compressibility with 
increasing guide field. The final particle energy spectrum in 
the largest simulation is best fit by the inclusion of new thermal 
components at temperatures 2.1 m e c 2 , and 3.5 m e c 2 , in addition 
to the initial thermal component with temperature m e c 2 . We, 
however, acknowledge that a larger size or longer duration sim- 
ulation is likely to produce a still more pronounced energized 
component, possibly even a population described with a power 
law spectrum. 

Energetic positrons (electrons) with Lorentz factors 7 > 
30 are moderately beamed in (opposite to) the direction of 
the initial current flow with median inclinations of ~ 30° - 
40°. The degree of beaming is determined by a particle's 
energy gain during acceleration. We speculate that more highly 
magnetized plasmas and reconnection sites with larger size 
X-line regions should give rise to stronger beaming. 

In this work, we have investigated a narrow range of magne- 
tizations with a ~ 0(1), but astrophysical reconnection sites 
can also have high magnetizations a 3> 1. We can speculate 
about the applicability of our results in the latter limit. The 
linear tearing mode responsible for the initiation of reconnec- 
tion is insensitive to the degree of magnetization far from the 
current sheet. The phase decoherence that produces the initial 
breakdown of translational invariance is determined by the 
tearing mode growth time and should thus also persist at high 
magnetizations. Therefore, we expect the qualitative structure 
of the reconnection region at higher values of a to be similar 
to that found in our simulations. 

The primary effect of high magnetization is that the Alfven 
velocity approaches the speed of light, which could give rise to 
ultrarelativistic outflows from the X-line region. In such out- 
flows the inertial term of the generalized Ohm's la w becomes 
impor tant in the breaking of flux freezing ( |Hesse & Zenitani 



2007). This in turn may increase the dimensionless reconnec 
tion rate r rec relative to the value found in our simulations. It 
remains to be seen whether the associated reconnection pro- 
cess is more or less intermittent. An increased magnetization 
is likely to increase the efficiency and the degree of beaming 
in particle energization. 

We consider these results and the immediate questions they 
raise an incremental step in the development of a multiscale 
view of collisionless plasma self-organization during magnetic 
reconnection. Further work is clearly required to place our 
key finding, the evolution of the simulated, periodic reconnec- 
tion layer into a disordered network of interacting magnetic 
flux ropes, in the macroscopic context of a realistic reconnec- 
tion site characterized by outflow boundary conditions and 
altogether different field line asymptotics at large distances 
from the X-line. It will be particularly interesting to see if 
the reconnected-flux-carrying outflow from the macroscopic 
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X-line will possess the disordered, interlinked magnetic field 
topology we observe and investigate what will be the character 
of magnetic fluctuations in the outflow. 
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